function diff = eta_fun(m,param)

    % hiring rate eta(m) 

    if param.zeta==0
        diff = - param.sigma^2./2./(1-param.alpha).*m.*delta_p_fun(m,param)./delta_fun(m,param);
    elseif param.K>0
        % set hiring rate in replacement region
        diff=zeros(size(m))+(param.zeta+delta_fun(m,param)).*(m>=param.m_h).*(m<=param.m_e);
        if min(m)<=param.m_h
            % set hiring rate in natural wastage region
            diff((m<=param.m_h)) = zeros(size(m(m<=param.m_h)));
        end  

        if max(m)>param.m_e
            % set hiring rate in expansion region
            diff((m>param.m_e)) = - delta_p_fun(m(m>param.m_e),param).*q_fun(m(m>param.m_e),param)./q_p_fun(m(m>param.m_e),param);
        end    

    else
        diff=zeros(size(m));
        if max(m)>=param.m_h
            diff((m>=param.m_h)) = - delta_p_fun(m(m>=param.m_h),param).*q_fun(m(m>=param.m_h),param)./q_p_fun(m(m>=param.m_h),param);    
        end  
    end

end